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ABSTRACT 

The term 'trajectory problem' is taken to include problems that can arise, 
for instance, in connection with contour plotting, or in the application 
of continuation methods, or during phase-plane analysis. Geometrical 
techniques are used to construct difference methods for these problems to 
produce in turn explicit and implicit cii ;ularly exact formulae. Based 
on these formulae, a predictor-corrector method is derived which, when 
compared with a closely related standard method, shows improved performance. 
It is found that this latter method produces spurious limit cycles, and 
this behaviour is partly analyzed. Finally, a simple variable-step 
algorithm is constructed and tested. 


•Visiting Scientists at the National Research Institute for Mathematical 
Sciences, CSIR, P 0 Box 395, Pretoria, South Africa. 
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1. INTRODUCTION 

We consider an initial-value problem for the autonomous system of ordinary 
differential equations 

dy 

3f*f(y). (1.1) 

where y is a vector in , 

In a number of practical applications the interest lies in obtaining the 
curve traced by the solution y(*) rather than in finding the actual cor= 
respondence between values of the independent variable or parameter t 
and points on that curve. These applications include the computation of 
trajectories in mechanical problems, the plotting of the phase-plane of 
second-order autonomous differential equations [2], and the study of 
solution fields of nonlinear equations 11,51. We shall employ the terra 
trajectory problem to refer to these cases. 

By definition a trajectory problem is not altered if the independent 
variable in (1.1) is replaced by a new variable u = irlt), where v is 
differentiable and monotonic. On the other hand the performance of a 
numerical method when applied to (1.1) depends heavily on the particular 
parametrization [61. To overcome the difficulties associated with the 
choice of this independent variable, the following devices come easily 
to mind. 

(i) Use of one of the coordinates, the first say, of y as independent 
variable. This procedure reduces by one the dimension of the sys= 
tern, but suffers from the disadvantage that the integration cat not 
be carried beyond a point y for which f^(y) = 0. It should also 
be noted that this procedure is not invariant with respect to 
rotation of the axes in the y-space. 
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(ii) Parametrization of the curve by its unique intrinsic parameter, i.e. 
its arc length s. This is equivalent to replacing (1.1) by 

| = == 

since now Idy/dsI = 1. (Me shall here not be concerned with singular 
points where f(y) * 0.) The use of the arc length and some of its 
modifications has been advocated by H B Keller L4] in the context 
of the solution of nonlinear equations. See also [6]. 

For the two-dimensional case (m=2) Lambert and McLeod [2] have intro* 
duced a sucressful modification of the idea in (i). They use the mid= 
point rule by rotating locally the axes in the y plane so as to have the 
tangent to the solution at the latest computed point playing the role 
of positive direction of the independent variable. This local rotation 
renders their method intrinsic in the sense that it does not depend upon 
the orientation of the axes in the y plane. Lambert and McLeod prove 
their mettiod co be circularly exaci , i.e. if the trajectory is a circle 
all the computed points will lie on the circle, provided that the starts 
ing points do and that no round-off error is present. Laurie [3] has 
extended the idea of local rotation to higher-dimensional equations. 

It appears to be desirable that a method should be circularly exact, as 
any m-dimensional curve can be approximated to second-order terms by its 
local circle (see Section 3). 

This paper continues the study of difference schemes specifically derived 
for trajectory problems. 

In Section 2 we present a simple geometrical way of constructing such 
methods. 

The local accuracy of the schemes is investigated in Section 3. 
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In Section 4 we define a circularly exact, fixed-step predictor-corrector 
algorithm that is closely related to the standard predictor-corrector 
method comprising the mid-point and trapezoidal rules in PECE mode. 

When both algorithms are tested in a number of problems the standard 
method is found to produce spurious limit cycles in some cases. It is 
proved that for a model problem the spurious cycles are local attractors. 

In the final Section we present a variable-step version of the circularly 
exact algorithm, a version whose step control strategy is based on a 
Milne device. Numerical examples are given. 

2. GEOMETRICAL CONSTRUCTION 

We illustrate the general idea by constructing the circularly exact method 

of Lambert and McLeod. This is an explicit, two-step formula which com= 

piites terms of the back points y^, y^^j and the back slope 

f - ^(y i)- We note that the points y , y , and the vector f_ , 

-n+1 -'in+1' tn Tn+l ,n*l 

uniquely determine a circle C^ in the m-dimensiunal space, the circle 
degenerating to a straight line if y^^j - y^, f^^^ are parallel. Choice 
of any point y^^^ makes the formula circularly exact. In particu= 

lar we can define y^^^ point on C^ sucli that ■ 

*~n+l ' ^n* depicts the two-dimensional plane 

spanned by the points y^. y^^j and the vector f^^j)- 

It is clear that with this choice 

-V2 ^ 'f'nrlLVl ' .rn’l ^+1’ 

where = ^n+1'^* ^n+1* ’ precisely the Lainbert-McLeod method 

as written by Laurie 13 1. 

3y construction the method generates points such that *yp^; - y^' is 
constant. In a 'variable-step' implementation one may wish to increase 
or decrease the Euclidean distance between consecutive points, and this 
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can be achieved by changing the choice of on C^, as will be shown 
in Section 5. 

Turning now to the general idea, suppose that we are given a fwily of 
curves such that an individual mender of the family can be determined by 
M linear conditions (when m • 2 three conditions determine a circle, 

four a parabola, five a general conic, etc ). Then M pieces of in» 

formation from the back data can be used to determine a cur^e of the 
family, and any choice of the next point on this curve will yield an 
explicit method which is exact whenever the trajectory belongs to the 
given family. 

This idea can also be employed to derive implicit methods. In this case 
the slope at the next point appears in the formula, and only M - 1 pieces 
of information from the back data are required. As an illustration, let 
us derive a circularly exact one-step method. From Figure 2 we see that 
when the solution is a circle, y^^j - y^ bisects the angle between the 
unit vectors F^, F^^j. 

Therefore 

>n+l ■ ^n ' 

where k is a parameter, yields the method sought for. Of course (2.2) 
is nothing but the trapezoidal rule applied to (1.2) with step-s^ze k. 

3. THE TRUNCATlOh ERRO R 

In this section we attempt to define the concept 'truncation error' for 
methods such as (2.1). In order to motivate the definition, let us con- 
sider first the formula (2.2). When this is viewed as the usual trape« 
zoidal rule applied t. (1.2), the standard procedure is to define the 
truncation error at a point y(s^) of the trajectory by 
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A Taylor txpansion revtals that as k ■» 0 

TE - - L k" ^ ♦ 0{k“) (3.2) 

I* *o 

and accordingly one says that the method Is of second order. We recall 
that If we denote by t, n, b the local tangent, (first) normal and second 
normal unit vectors, respectively, the derivatives of y w.r.t. s can be 
expressed as follows: 

y . t (3.3) 

• • 

y ■ icn 

• • • • j 

y • icn - If t ♦ Kib. 

Here a dot represents differentiation with respect to the arc length s 
and K, T the first and second curvatures. When the curve Is three-d1< 
menslonal the terms binomial and torsion are often used to refer to b 
and T respectively. 

From these expressions we see that In the neighbourhood of a point any 
m-dimenslonal curve can be approximated to second-order accuracy by the 
circle which shares Its curvature, and tangent and normal vectors. When 
(3.3) Is taken Into account (3.1), (3.2) can be written as 

TE . y(s^*k) - y(s^) - | [t(s„) ♦ t(s^>k)l . (3.4) 

* * T? [/(Sq) n(s^) - <^(s^) t(sj ♦ •^(So)"(So)b(s^)l ♦ 0(k'*). 

When the true trajectory Is a circle, i(,t e 0 and (3.4) becomes 

TE • kVt ♦ 0(k“). 





(3.5) 



Th# f»ct th*t we ere eleeling with a circularly exact »ethod 1* not appa* 
rent fro» (3.4). This is due to th# fact that th# truncation error 
locally Rieasures th# distance between the conputed point y^^^j and th# 
exact y(s^^j) (when y^^ • >{*„)). whilst we are Interested in the distance 
between and the trajectory. 

As an alternative we shall define the concept of reduced truncation error 
(RTE) which has the followinj property: whenever the method is exact 
for a family of curves in the sense of the previous section, the RTE for 
a trajectory on that family vanishes identically. 

For the particular case of the trapeioidal rule we proceed as follows: 
we denote by h « h(k) the Euclidean distance between y^^j and y^ when 
y « y(s„), and then define the RTE at y(s„) by 

RTE . y* - y(s^) - ^ Lt(s^) ♦ t»l (3.6) 

where y* is the point on the trajectory such that ly’ - y(Sp)> • h and 
t* is the unit tangent vector at y*. 

Thus whenever a Step of the trapezoidal rule starting from y^ ■ y(s^) 
leads to a point y^^j which lies on the trajectory, we shall have y* • y^^^j 
and het.ce RTE • 0. 

Let US now expand the RTE (3.6) in powers of h. In order to do so we 
reparamrtrize the trajectory in the neighbourhood of y(s^^), taking as 
new parameter the Euclidean distance h($) • ly($) - Taylcr ex* 

pension of y(s) - y{i^) and use of (3.3) reveal that 

^ ^ ( 3 . 7 ) 

Now the standard rules for the differentiation of inverse and co^site 
functions yield the following expressions for the derivatives of y w.r.t. 


s 

dy/dh - t, (3.8) 

d^y/dh^ ^ trn, 

dVdh^ « in - 3/4 ♦ «:t b. 

Analogously, for the dcrlvttlvos of t one hes 

dt/dh • rn, (3.9) 

d^t/dh'* • ten - <c't ♦ iftb. 

Next, we eliminate k from (3.6), noting that 


Z'y' ' y(So)' 


2h 


(3.10) 


Substituting (3.10) into (3.6) and expressing the result in terms of the 
parameter h, we have 


RTE . y(h) - y(0) - ft(0) ♦ t(h)1. (3.11) 


We now Taylor-expand, using (3.8), (3.9) to replace the derivatives of 
y, t, and arrive at 


rTE • - ~ h’frn ♦ nb1 ♦ 0(h“). (3.12) 

We note that the curvature does not appear alone in the leading terms 
of the RTE, in ag re e m e n t with the fact that RTE « 0 if r, i ? 0. (In 
fact it can be shown that the whole Taylor series for RTE does not in- 
volve terms which contain only the curvature.) 


The idea we have just illustrated in the case of the trapezoidal rule 
can be extended to ocher menders of the class of methods introduced in 
the previous section. For instance for the method (2.1) one would de- 
fine 


RTE • y»’ - y(s„) ♦ 2 1 t*^ (y‘ - y(s^))l t* . 


(3.13) 
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i*h#re y*, y” »rt th# points on th# curve such titet •y(*Q)-y*< • ly*-y**l » h 
with h equal to the CMtttant dittance between any two consecutive points. 

He now find 

RTE ■ ^ ten ♦ ttb] ♦ 0(h'*). {3. 14) 

This Idea of an RT£ can be employed to derive estimates of the global 
accuracy of the methods. The details will be given elsewhere. 

4. A C 1RCUUUU.T EXACT PREDICTOR-CWRCCTW tffTIffl) 

CQeq>ar1son of (3.12) with (3.14) shows that the Implicit circularly 
exact metiwd (2.2) has a smaller error '.onstant than the explicit method 
(2.1). Therefore It Is reasonable to consider the Idea of combining 
the two methods In a pred1ctor*corrector pair. He suggest the following 
formulae: 

?ihZ ‘ Untl * I 1 

"n+1 . n+2' 

where h - ly^ - >jl . 

Note that 'y ^^2 - y^|' » »nd that the step-length k of the 

corrector (2.2) Is changed from one step to the next In order to guaran< 
tee that ' J!nel' * »’• 

Hhen the trajectory 1$ a circle and y^^, y„^j. He o" the trajectory, the 
predictor yields a point y^„^j on the circle with ly^^j * y**n^ 2 * * *>• 
Therefore y^^j • *nd the method Is circularly exact. 

Formulae (4.1) were tested in several numerical examples, and In order 
to establish a fair coa^>ar1sor the following method was used; 
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?’m! ■ ?n • ?».!• " ’I 

• '«1 • l*«l (?■»! • f’r.z'- 

l.e. th€ pr*d1ctor-corr«ctor Method bated on the Mid-point end trepetol* 
dal rules used In PECE Mode. Recall that F * f/lfl. 

It should be stret: .d that if correction to conver 9 ence rather than the 
PECE Mode had been used, one kould have had the circularly exact Method 
(2.2). Itowever (4.2) Is not circularly exact, at will be clear froM the 
following discussion. 

Suisse that (4.2) It applied to the two-dlMenslonal problcM 

fj • • >2 (*-3) 

^2 ■ >1 

whose trajectories are circles centered at the origin. 

This problei) 1* best analyzed by means of polar coordinates. NaMely let 
us describe each of the vectors y^ generated by (4.2) by the radius 
p • ly I and the angle formed by y_ ,,y_. Then, after some tsanipu' 

f| Ifl n -n*i -n 

latlon, it Is found that y^^^ f'. obtained from by ineans of the 

forMulae 

<^n*2 * ’ 

cos - (P „„ ♦ p -„,2 - k- cos'B)/(2 p„^ 2 ). (4.4b) 

where B Is a function of k, o^, given by 

cot 2ti • p^ cos a^^j/(2k - si na^^j). (4.5) 

We see from (4.4a) that in general the radius r does not reaain constant 
for all iterants and therefore that the method is not circularly exact. 

It is useful to take this discussion further as follows. Formulae (4.4) 


describe a tMO-step recurrence for the cosfMJtatlon of 
teres of «n+j) *"<* P^' possible to reforwiate this recur* 

rence as one having only one step, by Increasing the diaenslon of the 
vectors Involved. Naaely with r^^ ° arrive at the recurrence 

W = Pn4l* 

‘^n+2 ° 

®n+2 ° ^P^n+1 ' •‘^cos^e)/(2 P„+i «)]. 

where 6, R satisfy 

cot 26 = r„^j cos a„^i/(2k - r^^j sin (4.7) 

R * (k^ cos^S + - k sin 26) ^ 

Now (4.6) describes the transformation of (r^j^j. P^j» “^+1^ 
^'^n+2’*’n+2’“n+2^‘ easily verified that (k/2, k/i.-n/Z) is a fixed 

point of this iteration. 

We conclude that if (4.2) is applied to the model system (4.3) with 
lygl = *>j* = l'/2 and yjj, y^ forming an angle of ir/2, then each subse* 
quent iterant also lies on a circle of radius k/2 and is tt/ 2 radians 
from the previojs iterant. We shall use the term 'spurious limit circle' 
to refer to this circle of radius k/2. 

The Jacobian matrix of the transformation (4.6) evaluated at the fixed 
point is found to be 

0 1 0 

0 0 -k/6 

0 -2/k 0 

with eigenvalues 0, - /3/3. Since these are srnaller in magnitude than 
unity, the fixed point is a local attractor, i.e. initial vectors Yq.V, 
near the spurious limit circle and forming an angle near to ir/2 will pro* 
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d'ice a sequence of iterants which converges to the spurious limit circle. 
In fact we shall see In what follows that the Iterants converge to that 
circle even If the initial vectors are far from it. 

We are now in a position to report several m«erical tests on methods 
(4.1), (4.2). In all the examples the 'exact' trajectory was calculated 
employing the usual fourth-order, fourth-stage Runge-Kutta method, which 
also provided the additional starting value. 

When using the Runge-Kutta method, a step-size one-tenth that of the 
predictor-corrector algorithms was taken. In the figures the points 
produced by the Runge-Kutta method have been joined by a continuous curve, 
those produced by (4.1) being indicated by circles 'O' and those produ= 
ced by the trapezoidal rule indicated by crosses 'X', and joined by a 
broken line for additional clarity. 

As a first example we consider the model problem (4.3). The initial 
point was (0,1) and the step 1. ('Step' means, of course, h in formulae 
(4.1), k in formulae (4.2).) The results have been plotted in Fig. 3. 
Ninety-eight points were computed for each algorithm. Those correspond^ 
ing to the circularly exact method fall repeatedly on the inscribed hexa= 
gon, showing numerical stability. The points corresponding to (4.2) 
spiral very rapidly towards the spurious limit circle, and from the six= 
teenth onwards lie on that circle (within the accuracy of the plot). 

Fig. 4 corresponds to the same problem and initial condition, but the 
step is now 0.37. Note that the radius of the spurious limit circle 
has decreased, in agreement with our earlier discussion. 

The second example is the system 

= - i'2. (4.8) 


f^ = sin y^. 
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which is equivalent to the well-known pendulum equation. 

The initial point was (0,1) and the step 0.5. The behaviour of the me= 
thods was very similar to the one we have seen in the first example. The 
points produced by the circularly exact method were, within the accuracy 
of the figure, on the exact integral curve. The solution given by the 
method (4.2) spiralled in and reached a limit circle of radius 0.25, far 
from the true orbit. This value of the radius is precisely that of the 
spurious limit circle for the model problem. This is no surprise as the 
phase-planes of (4.3), (4.8) near the origin are very similar. 

The next example is the van der Pol system 

fl = /2 - -l(yi - 3yi) (4.9) 



The results illustrated by Figs. 5 and 6 both refer to a step 1.5 but the 
starting point was (10,10) for the former and (0,1) for the latter. We 
see that in both instances the circularly exact method identifies correct* 
ly the limit cycle of the system, whereas the results given by the 
method (4.2) suggest a 'spurious' limit cycle whose diameter is roughly 
half the true one. Neither method does well in the descending section 
of the trajectory in Fig. 5. We shall see later that the integration of 
(4.9) is comparatively difficult in that region. 

For Figs. 7 and 8 the step was 1. Again the method (4.2) produces a 
spurious limit cycle. It appears that the size of the spurious limit 
cycles obtained does not depend on the initial point but only on the 
step size. 

The last example had 

r’l « yg (2 y^ + y^) (4.10) 
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and initial point (0,1). The results for h . k - 0.5 are depicted in 
Fig. 9. The points corresponding to (4.1) are reasonably close to the 
true trajectory even when five orbits have been completed, while the 
method (4.2) once more yields an incorrect picture of the situation. 

We conclude that for the problems considered the geometrically derived, 
circularly exact algorithm (4.1) is better suited than its standard 
counterpart. 

5. VARIABLE STEP 

In this Section we construct and test a variable-step version of the 
circularly exact method (4.1). It should be emphasized that our aim is 
to demonstrate the possibility of such a construction rather than to 
develop a sophisticated code. 

We first derive a variable-step circularly exact predictor formul*. 

Given y^, y^^j* and a positive number h^^j, this formula will yield 
the point satisfies •y^„^2 ' ^n+1* * ^n+1 

circle determined by y^, y^^j, F^^j- Fig. 10 depicts the two-dimen= 
sional plane defined by the points y^, y^^^ and the vector F^^^. If we 
denote by y the angle between y^^j - y^ and F^^j, then the central angle 

P 

subtended in by y^, y^^^j is 2y. Therefore the angle between y,, * V n+2 
and y^^j - y*’p^2 ^^at an inscribed angle is equal to one 

half of the corresponding central angle.) 

Next let 6 be the angle between F^^j and y^^^2 ' ^n+l" angle 

between y^^j - y^ and y^„.^2 ‘ ^n+1 is and con, deration of i.he 

triangle with vertices y^,. y ^^2 conclusion that the 

angle between y^^j • y^ and y’’p ,^2 ■ y„ is also 6. We have denoted by 
N^^j the unit normal vector to at y^,,j- 

We are now in a position to derive the required formula. We project 

y"n.2 - ^n.l ^n+1- ?n.l 
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■ >n+l “ ^+1 « !n+l ^ ^+1 « ?n+r 

The 6ram-Schm1dt procedure enables us to express the normal vector !!n.l 
in terras of *s follows: 

Nn+1 ■ > ^n^l * ('’n 

where h = ly , - y I. Next 6 can be eliminated by use of the sine theo= 
n tn+i tn 

rem in the triangle y^, = 

hj^/sin Y * 

Finally y is related to F^^j, y„^j. y^ by the formula 

fn*l {^041 ■ yn> ' ^ 

When (5.2), (5.3), (5.4) are substituted into (5.1) the following predic* 
tor formula is obtained: 

y'*n+2 * ^n+1 ^ !n+l ^ U ‘ ^ntl^ {^-5) 

where 

O . h“ , 

(5.6a) 

Bn = !Vl(>n.l-?n)- <5.6b) 

Formula (5.5) reduces to formula (2.1) if h^^j « h^. It should also be 

noted that y **„^2 *’* defined ^.1 chosen larger than the 

diameter d of C . From Fig. 10 this diameter is h /sin y, whence using 
n n ^ n 

(5.4), (5.6b) we obtain 

d„ . hj /(»; - BJ*. (5.7) 

In fact the algorithm we shall describe later imposes the condition 


The corrector formula is written in the form 
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t='l 

so thot ly„j - y„„l . h„„. 

In order to control the step-size a Milne device can be employed. Let 
y* be the point such that ly* - y^^jl • h^^j and y* lies on the trajectory 
through Then as in Section 3 

y* - y'’n.2 = . h-;^j h^)[^n . KTb], (5.9) 

y* - w " 

and elimination of the term in square brackets leads to 

y* - yn.2 = Vi * 2 ^)]fyn.2 ' /n.2l' (^.ll) 

We considered the following algorithm 

(1) Given yo^yj’^O’^l' ^ ^ ^0 ' 'yi ’ yQ* * 0; 

(2) Evaluate F^^j. Use (5,6b), (5. 8) to compute B^, 1/d^. If l/hj^^j< 2/d^, 

(3) Compute y*^^^2 to (5.5), evaluate and form y^^^ 

(formula (5.8)). 

(4) Use (5.11) to estimate the error e = ly* - yp+2* • ^n+1 ' ^n+1 

(e/e) j; 

(5) If e > E , set h^^j ' h*^j and go to (3); 

(6) Print y^^^' '’n+2 * *’n+r " ' 5° to (2). 

The trapezoidal rule in correction-to-convergence mode was used to com= 
pute yj and initialize the algorithm, which was tested with several 
tolerances e and various initial points in the systems (4.9) , (4. 10) . 


The following three-dimensional system was also considered; 
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fj • yj. (5. 12) 

h “ ^ 1 * 

f3 « 4 yj y^. 

Fig. 11 shows the results for tne system (4.9) with e ■ 0.001 and 
yg * (30,30). The true trajectory starting from (30,20) is also dep1c= 
ted in order to display the rapid convergence of the integral curves in 
the vicinity of the vertical portio'. C,D. It is well-known that this 
convergence forces any explicit algoritNn to take a small step. By com> 
parison the step is larger along AB, where the neighbouring integral 
curves are almost parallel. 

Fig. 12 also refers to the system (4.9), but now e = 0.005 and y^ • (0,1). 
The maxinHjm Euclidean distance between consecutive points is 1.6. 

Fig. 13 corresponds to the system (4.10) with e • 0.001 and y^ * (0,1). 

The system (5.12) was integrated starting from (1,0,1). ■ The true solu« 
tion is given in parametric form by 

y^(t) * cos t (5.13) 

=-sin t 

yj(t) * cos 2t 

The integration was stopped when roughly a quarter of an orbit had been 
completed. This corresponds to an arclength of 2.63. The curvature is 
initially 2.0, decreasing to 0.1 and increasing again to 2.0. When the 
toleranca was 0.01, nine steps were taken and the final point lay at a 
distance ot 0.025 from the true integral curve. When the tolerance was 
decreased to 0.0001, thirty-two steps were required and the final error 
was 0.003. 


We wish to emphasize that the algorithm presented here can be easily adap 
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ted to yield several geometrical elements of the trajectory such as 
tangent and normal vectors, curvature, arc length, etc. 
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